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ABSTRACT 

We have performed three-dimensional, non-LTE (non-Local Thermodynamic Equilibrium) radiative 
transfer calculations for 12 CO and 13 CO lines, applying them to our high-resolution hydrodynamic 
models of the 'torus' around a supermassive black hole in an active galactic nucleus. The hydrodynamic 
simulations reveal inhomogeneous and turbulent gas structure on a sub-pc scale in a circum-nuclear 
starburst region. Thick disks interlaced with filaments, clumps and holes are naturally formed due 
to the interplay among energy feedback from supernovae, self-gravity of the gas, galactic rotation, 
and radiative cooling. The intensity maps of the molecular lines for the circum-nuclear disks show 
a clumpy structure reflecting the intrinsic inhomogeneity and turbulent motion of the gas disk. The 
fine structure of the 'torus' could be resolved in the nearby active galaxies using the Atacama Large 
Millimeter/submillimeter Array. We also found that the CO-to-H2 conversion factor (X-factor) is 
not uniformly distributed in the central 100 pc region. The X-factor derived for 12 CO (J =1-0) 
intensity depends strongly on the intensity, whereas the X for 12 CO (J =3-2) is nearly constant 
over two orders of magnitude of the intensity. The suggested conversion factor for the molecular gas 
mass is A"co(j=3-2) ~ 0.27 x 10 20 cm -2 (K km s^ 1 ) -1 . The line-ratio of high-J transitions of CO 
is not uniformly distributed in the nuclear disk, and the apparent ratio depends on the beam sizes. 
'Observed' 12 CO (J=4-3)/ 12 CO (J=2-l) can differ from the intrinsic ratio by maximally 20%. 
Subject headings: ISM: structure, kinematics and dynamics — galaxies: nuclei, starburst — method: 
numerical 



1. INTRODUCTION 

The galactic centers of spiral galaxies are important 
regions for the evolution of galaxies. It has been widely 
believed that most 'large' galaxies harbor supermassive 
black holes (SMBHs) at their centers, and also that the 
observed scaling relation between masses of black holes 
and the galactic bulges suggests that they have evolved 
toget her during the hierarchical f ormation of galaxies 
(e.g. iKauffmann fc Haehneltl 120001) . The galactic cen- 
ters are often very active. Nuclear starbursts and vari- 
ous kinds of active galactic nuclei (AGNs) are observed 
in nearby galaxies and also in galaxies at high redshifts. 
These active phenomena must be related with the evo- 
lution of the SMBHs and galaxies, as well as the inter- 
galactic medium through galactic winds and radiative 
feedback, especially at high redshifts. All these phenom- 
ena are related to the structure and dynamics of the in- 
terstellar medium (ISM) in the central 100 pc regions 
of galaxies. For example, AGN activities are supported 
by gas fueling, whose mechanisms remain unclear, espe- 
cially in the inner 100 pc from the central engine (see 
Wada 2004, and references therein). We have not yet 
determined what triggers the nuclear starbursts that are 
conn ected to the formation of st ellar cores and the SMBH 
(e.g. iNorman &: Scovillel Il988j) . However, in the next 
decade, we may expect drastic change in our observa- 
tional knowledge of the ISM in the circum-nuclear region 
of galaxies. The ~ 0.01" spatial resolution with the next 
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generation millimeter and submillimeter interferometer, 
ALMA (the Atacama Large Millimeter Array) may re- 
veal a sub-pc structure of the ISM. We now need to im- 
prove our theoretical modeling of the ISM in the galactic 
central regions in order to make direct comparison with 
future observations. 

The ISM in galaxies is inhomogeneous and multi- 
phase, and this is especially important in studying the 
gas dynamics in circum-nuclear regions, because the typ- 
ical scale of the inhomogeneity of the ISM (~ 10 pc) is 
not negligible compared to the size of the system (~ 100 
pc). The isothermal approximation (e.g. the sound ve- 
locity ~ 10 km s _1 ), which has often been used for mod- 
eling gas dynamics on a galactic scale, is not adequate 
for the nuclear gas dynamics, particularly when energy 
and radiative feedback from stars is taken into account. 
The multi-phase nature is also important for comparing 
numerical models with observations. Radiative trans- 
fer calculations, for continuum and various lines emitted 
from the AGN regions, require realistic density, temper- 
ature and velocity distributions of the ISM as input. In 
order to achieve realistic numerical models of the ISM, 
sub-pc spatial resolution for diffuse gases as well as for 
dense gases and resolving strong shocks is also essential. 
We should also take into account the self-gravity and 
three-dimensional dynamics of the ISM. 

Instead of adopting phenomenological approaches for 
the multi-phase ISM, in which man y free parameters 
and assumptions are unavoidable (e .g. iCombes &: Gerinl 
ll985j)lWada fc Normanl pOOl 12002) obtain global, inho- 
mogeneous models of the ISM with a sub-pc scale resolu- 
tion, by solving the governing basic equations of the ISM 
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with a high-accuracy Euler mesh code. In the present pa- 
per, using these results (three-dimensional density, tem- 
perature and velocity field data), we have made the first 
attempt to derive the molecular line intensities emitted 
from the nuclear starburst region with a supermassive 
black hole applying our non-LTE (non Local Thermo- 
dynamic Equilibrium), fully three-dimensional radiative 
transfer calculations. 

These kinds of 'synthesized observations' based on 
radiative transfer calculations has rece ntly been ap- 
plied to turbule nt mo lecular clouds (iPadoan et al.1 

199St IQssenkonfl B IQssenkonf fc Mac Low! 12002^ 
Pado an et alJ l|1998fl performed non-LTE radiative trans- 
fer calculations of C O and CS line s using a Monte-Carlo 
integration scheme l|Juvelalfl997(l for numerical results 
of the magneto-hydrodynamic simulations of turbulent 
clouds. They found that the synthesized maps and spec- 
tra bear resemblance to re al mo lecular clouds in many 
respec ts. lOssenkopl l|2002|) and IQssenkopf fc Mac Low! 
(2002) made similar analysis of turbulent clouds, but 
they assumed a local LVG (Large Velocity Gradient) ap- 
proximation and isotropic radiation field. Ossenkopft 
et al. claimed that their simplified analysis is accu- 
rate enough for comparison between the observations of 
molecular clouds and the numerical models of turbulent 
clouds. These studies focused on the molecular clouds in 
disks of spiral galaxies, like our Galaxy, therefore we can- 
not simply apply their results to the interstellar medium 
in the nuclear starburst with a supermassive black hole. 
Apparently we should not assume isothermal, isotropic 
turbulence in the galactic nuclear region. Since the local 
turbule nt motion is of the sam e order as the global ro- 
tation l|Wada fc Normanll2002fi . a fully consistent treat- 
ment of the radiative transfer is necessary. 

In §2, we describe the numerical method and models. 
In §3, the typical structure of the ISM in the circum- 
nuclear region obtained by hydrodynamic simulations is 
briefly described. Then, 12 CO and 13 CO line intensity 
calculations are shown. The CO-to-H2 conversion fac- 
tor, the so-called "X" factor and the line-ratios are also 
discussed. A brief summary and discussion are presented 
in §4. 

2. NUMERICAL METHOD AND MODELS 

2.1. Hydrodynamic models 

The numerical methods for obtaining three- 
dimensional density, temperature and velocity fields of 
the gas around a SMBH, for which the line-transfer and 
excitation of mole cules are calculated, are the same as 
those described in IWada fc Nor man (2002). Details of 
the two-dime nsional version of the n umerical scheme is 
described in IWada fc _Ngnngrj_j20 i 0J^. In addition to 
the model described in IWada fc Normanl l)2002lk we also 
performed a model with a moderate supernova rate and 
a less massive black hole (see below). Here, we briefly 
summarize them. We solve mass, momentum, and 
energy conservation equations with the Poisson equation 
numerically in 3-D to simulate the evolution of a rotating 
ISM in a fixed gravitational potential. Here the external 
potential force in the momentum conservation equation 
is V$ ox t + V$bh i where the time-independent external 
potential is $ cxt = ~(27/4) 1 / 2 v 2 /(r 2 + a 2 ) 1 / 2 with a 
core radius of a = 10 pc and a maximum rotational 



velocity of v c = 100 km s 1 . The central BH potential 
is $bh = -GM BH /{r 2 + b 2 ) 1 ' 2 with b = 1 pc. The 
mass of the BH is assumed to 10 8 M Q (model A) and 
10 r Mn (model B). Model A is the same model studied 
in IWada fc Normanl <|2002f) . We also assume a cooling 
function A(T g ) (5K <T g < 10 8 K) with Solar metallicity 
and heating due to photoelectric heating, Tuv and 
due to energy feedback from SNe, Tsn- We assume a 
uniform UV radiation field, which is ten times larger 
than the local UV field. We assume an equation of state 
for an ideal gas with a ratio of specific heat 7 = 5/3. 

The hydrodynamic part of the basic equations is solved 
by AUSM (Advection Upstream Splitting Method) (Liou 
& Steffen 1993). We use 256 2 x 128 (model A) or 256 3 
(model B) Cartesian grid points covering a 64 2 x 32 
pc 3 (model A) or 64 3 pc 3 (model B) region around the 
galactic center (the spatial resolution is 0.25 pc). The 
Poisson equation is solved using the fast Fourier trans- 
form and the convolution method. The initial condi- 
tion is an axisymmetric and rotationally supported thin 
disk with a uniform radial density profile (the scale 
height in z-dircction is 2.5 pc) and a total gas mass of 
M g = 5 x 1O 7 M (Model A) and M g = 1O 7 M (Model 
B). Random density and temperature fluctuations, which 
are less than 1 % of the unperturbed values, are added 
to the initial disk. 

Supernova (SN) explosions are assumed to occur at 
random positions on the disk plane. The average SN rate 
is ~ 0.8 yr" 1 (model A) and 0.08 yr" 1 (model B). The 
energy of 10 51 ergs is instantaneously injected for each 
supernova into a single cell, as a form of thermal energy. 
The three-dimensional evolution of blast waves driven 
by the SNe in an inhomogeneous and non-stationary 
medium with a global rotation is followed explicitly, tak- 
ing into account the radiative cooling. Therefore the evo- 
lution of the SN remnants, e.g. their duration and struc- 
ture, depends on the gas density distribution around the 
SNe, and we do not need to assume free parameters, such 
as the heating efficiency. 

2.2. Line radiative transfer calculations 

Next, we try to 'observe' the star-forming disk around 
a SMBH, calculated separately using the hydrodynamic 
method described in §2.1, with molecular lines. A 
snapshot of density, temperature and velocity fields of 
128 2 x 64 2 (model A) or 128 3 (model B) or cells [i.e. 
each cell is (0.5 pc) 3 ], averaged from the original hydro- 
dynamics simulations with 256 3 or 256 2 x 128, cells is 
used to obtain 12 CO and 13 CO intensities. We have per- 
formed three-dimensional, fully non-LTE calculations to 
obtain CO line intensity distribution. We use the den- 
sity, temperature and velocity data in a 256 3 data cube 
from a snapshot of a three-dimensional hydrodynamical 
calculation, and the accelerated Monte-Carlo method is 
applied. We have developed a three-dimensional code, 
based on Hogerheijde & van der Tak (2000), and this is 
optimized to our vector-parallel supercomputer. Here we 
briefly describe the method. The equation of radiative 
transfer is 

d ± = -I v + S v , (1) 

where I v is the intensity at frequency v along a partic- 
ular line of sight described ds, the source function S u 
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is defined as S v = jy/au with emission and absorption 
coefficients j v and a v . The optical depth is defined as 
dr = a u ds. The coefficients j v and a v for a spectral 
line are determined by absorption and emission between 
levels u and I with number density n u and nf. 

jf = ^n u A u ^{v\ 
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where <p(v) is a line-profile function peaked around the 
frequency is = (E u - E{)/h, and A u i,B tu , and B ul are 
the Einstein probability coefficients. The line profile <j){v) 
is given by 
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The level populations are determined by the equation of 
statistical equilibrium, i.e. 

ni[Y, k< iAi k + H k ^i(Bi k J v + Cik)} = (5) 

^kyin-kAki + J^k^ink{BkiJu + Cm), 

where J„ is the mean intensity of the radiation field, 
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(6) 



with v = \E k - Ei\/h. Equations |QJ, © and © 
are solved iteratively. The integral in eq. © is 
evaluated using the Monte Carlo approach. Instead 
of us ing the original Monte Carlo approach l)Bernesl 
1979), here we used the impleme ntation developed by 
Hoge rheiide. fc van der Takl l)2000|) . in which the integra- 
tion is performed along the rays incident on each cell from 
infinity. This approach has the advantage of convergence 
of the radiation field and level populations. 

For each grid cell, Nn ay rays are calculated. First, we 
choose the wavelength of photons v incoming to the cell 
at x according to equation JIJ . The central frequency of 
these photons is taken to be equal to vl [1 + n- v(x)/c], 
where vl, n, and v represent the central frequency of 
this transition in the laboratory frame, randomly chosen 
normal vector of the ray from x and the local velocity 
of the gas at x. For these photons, the absorption and 
emission due to the grid cell at x' on the path are calcu- 
lated using equations (2) and (3). In equation ijljl. we put 
V{\ = Vjj [1 + n • v(x')/c]_for the central frequency. Then, 
we integrate equation Q from the boundary, where the 
cosmic microwave background radiation is assumed, and 
obtain I„(x). After averaging the specific intensities of 
Njiay rays as equation © by adding respective intensi- 
ties, we finally solve the rate equation and obtain the 
level populations for the cell at x. This procedure con- 
tinues till the level populations of each cell converge. 

We solved ten excitation levels of 12 CO and 13 CO, 
and iVfiay = 200 for each cell is used. In most cases, 
the radiation field, the level populations, and the fi- 
nal intensity distribution of the lower-level (J < 6) 
molecular lines are converged in about 20 iterations. 
The collision coefficients for 12 CO and 13 CO are taken 
fromlGreen fc Thaddeusl dl97^ and[G reen fc Chapman! 
(1978). Ten energy levels of these molecules are included 
in the calculation. We have assumed constant abun- 
dances of 5 x 10~ 5 for 12 CO and 10~ 6 for 13 CO. The 
'micro-turbulence', a in eq.Q is assumed 1 km s _1 . 



Computational time for one typical run of these ra- 
diative transfer calculations is typically about 4 hours, 
using 32 processors of a vector-parallel supercom- 
puter Fujitsu VPP5000 (the peak performance is 9.6 
Gflops /processor) . 

3. NUMERICAL RESULTS 

3.1. Input models 

Recent high resolution, 3-D hydrodynamical simula- 
tions revealed that a typical three-dimensional density 
structure of the ISM in the central 100 pc is character- 
ized by a 'tan gled' network of filaments and dense clumps 
(Wada 2001). These filaments are a natural consequence 
of the gravitational and thermal instabilities, tidal in- 
teraction between dense clumps, and local and global 
shear motions. The probability distribution function of 
gas density in an equilibrium state is well represented 
by a single log- normal function (Wada & Norman 2001; 
Wada 2001). In other words, this reflects the fact that 
the denser media occupy smaller volumes, and the log- 
normal function implies that the inhomogeneous struc- 
ture is a consequence of highly non-linear ways. This 
nature of the density field would be related to the em- 
pirical relation between the global star for mation rate 
and t he average surface density of the ISM l)ElmegreerJ 
2002) . During the development of these inhomogeneous 
structures, a small part of the galactic rotational en- 
ergy is converted to maintain a steady turbulent motion 
of the ISM [see also Wada, Meurer, & Norman (2002) 
for two-dimensional cases]. The ISM around a SMBH 
(10 6-9 Mq) on a 10-100 pc scale is basically the same 
structure as mentioned above. However, the gravita- 
tional field of the massive black hole affects the scale 
height of the disk. Energy feedback from supernova ex- 
plosions, which is expected in the clumpy medium, also 
determines the global structure of the nuclear gas disk. 
Therefore, the scale height of the thick disk is determined 
from a balance between supernova heating and turbulent 
energy dissipation due to the radiative cooling; in a black 
hole and galactic gravitational field (jWada fc Normaiil 



Figure ^ is density distribution of model A (Mbh = 
10 8 M Q with a relatively high star formation rate ~ 
0.8M© yr _1 ), and it is clear that the gas around the cen- 
tral massive black hole is filamentary and clumpy, and 
the scale height is larger in the outer region. As pointed 
out in Wada & Norman (2002), the internal inhomoge- 
neous structure of the 'torus' is not steady, but rather 
seemingly turbulent, although it globally rotates around 
the center. The global flared disk shape does not change 
significantly during many rotational periods. The turbu- 
lent motion is enhanced by supernova explosions which 
are assumed to randomly occur in the disk, and their 
blast waves blow the gas up from the disk plane. 

Figure [5] shows the density structure of model B, in 
which a less active supernova rate (0.08 yr" 1 ) and a less 
massive black hole (10 7 M Q ) are assumed. The structure 
is similar to those seen in Figure ^ 

3.2. Molecular line-transfer calculations 

Since most parts of the gas mass in the torus are in cold 
(T < 100 K) and dense (n > 10 2 cm" 3 ) phases, we ex- 
pect that molecular lines, such as 12 CO (J =1-0, 2-1, ...), 



4 



Observing the obscuring torus 



are good tracers to observe the molecular torus. Here we 
demonstrate how the dense clumpy 'torus' is observedby 
the 12 CO lines, using three-dimensional radiative trans- 
fer calculations, without assuming the LTE described in 
§2. Figure shows channel maps of 12 CO ( J = 2 — 1) 
for the torus shown in Figure 1 with a viewing angle of 
45 degrees from the rotational axis. If we obtain channel 
maps like Figure 2 by future radio observations, we can 
confirm the ISM in the circum-nuclear region is globally 
rotating at ~ 200 km s -1 , and the internal velocity dis- 
persion for the cold and dense gas is as large as ~ 50 
km s _1 . The large random motion inside the torus sug- 
gests that the torus is geometrically thick, provided that 
the scale height of the torus is determined by a balance 
between vertical components of the centrifugal force due 
to t he central black hole an d vertical velocity dispersion 
(see iWada fc NorTnan1l2002]) . 

Figure 01 shows integrated intensity maps for model A 
in Figure H We convolved the maps with four different 
'beam' sizes to see how the difference affects the apparent 
structure inhomogeneity of the integrated maps. As is 
clearly seen in Figure 3b, the clumpy structure of the 
torus cannot be resolved with a beam size of about 1/10 
of the disk size, although the torus-like structure can be 
still observed. Suppose this torus is located in the Virgo 
cluster, 0.01 arcsec angular resolution, which is expected 
to be attained by the ALMA, approximately corresponds 
to the beam size in Figure 3c. Therefore, the clumpy 
structure of the molecular tori will easily be detected. On 
the other hand, the one arcsec beam, which is typically 
achieved by present-day radio interferometers, such as 
the Nobeyama Millimeter Array, is almost comparable 
to the torus size, therefore even the central hole of the 
torus could not be resolved, even if it in fact presents. 

The two panels of Figure are volume- weighted and 
mass- weighted histograms for a density and temperature 
in m odel B. As in the usual three-phases model of the 
ISM l|McKee fc Ostrikerll977|) . three stable branches can 
be seen in Figure 5a, i.e. cold, dense gases (10 < p < 
1O 5 M pc" 3 and T g < 10 2 K), warm gases at T g ~ 10 4 
K, and hot, diffuse gases at T g ~ 10 6 K and p < 0.1M© 
pc~ 3 . The temperatures of warm and cold gases are 
maintained by the photoelectric heating of the dust by 
the UV radiation, and shock heating in the compress- 
ible turbulence due to supernova explosions. They are 
roughly in pressure equilibrium, but one should note that 
the dispersion of the pressure is quite large for a given 
density or temperature, which suggests that structure of 
the multi-phase ISM may not be determined under a con- 
stant pressure. Figure [SJd, the mass- weighted phase dia- 
gram shows that the gas mass is dominated by the cold, 
dense phases. Mass fraction of the gas with T g < 50 K 
is about 20% of the total mass. 

Figure [5] shows line-ratios of 13 CO for the edge-on view 
of model B. The flared disk shown in Figure [3 is clearly 
seen with the ratio of 13 CO (J=1-0)/ 13 CO (J=2-l) - 2, 
and 13 CO (J=3-2)/ 13 CO (J=2-l) ~ 1, whereas in 12 CO 
line-ratios (Figure EJ, the flared disk is not clear. These 
results suggest that the 12 CO/ 13 CO line-ratio is not uni- 
form in the nuclear starbursts even for the constant abun- 
dance ratio that we assumed. This kind of fine structure 
of the molecular gas shoul d be detected in nearby Seyfert 
2 galaxies with starbursts <)Storchi-Bergmann et al.l200lJ: 



iCid Fernandes et al.l 12003^) by future observations with 
the ALMA. 

Figure |S1 shows distribution of the CO-to-H 2 conver- 
sion factor (X-factor) calculated from column density 
of the gas and integrated intensity of 12 CO (J =1-0). 
The starburst model (Fig. EJ) is observed from a view- 
ing angle of 45°, convolved with a Gaussian kernel of 
FWHM 2.5 pc (five grid points). It is clear that the 
X-factor is not uniformly distributed. In most regions, 
the X-factors are smaller than the value estimated for 
the local clouds in the Galaxy [2-3 xlO 20 cm" 2 ( K km 
s" 1 )" 1 ] l|Solomon. Rivolo. Barrett, fc Yahill fl987^ . Its 
distribution roughly follows the column density distribu- 
tion. This is also clear in Figure[7| in which the X-factors 
for four lines of 12 CO are plotted against the integrated 
intensity. The X-factor of 12 CO (J =1-0) depends on 
the intensity, and the larger X-factor (i.e. higher den- 
sity) is obtained for stronger intensities. In other words, 
the column density of the molecular gas in dense regions 
may be under-estimated if a constant conversion factor 
is assumed. On the other hand, the X-factor of 12 CO 
(J =3-2) is nearly constant for an integrated intensity 
of less than 2000 K km s _1 . These results suggest that 
12 CO (J =1-0) is not a good tracer for gas column den- 
sity and total gas mass for the dense nuclear gas disk 
with the starburst, but 12 CO (J =3-2) can be used to 
estimate the column density of the molecular gas in a 
wide range of intensity. The average X-factor for 12 CO 
(J =3-2) is (0.27 ± 0.04) x 10 20 c m - 2 (K km s' 1 )- 1 in 
the range of 10 2 <J co <2x 10 3 K km s" 1 . These 
results are consistent with the well known fact that the 
mid-J CO transitions are suitable for derivation of the 
column density of warm (e.g. 30-50 K) mol ecular clouds 
in our Galaxvfe.g. lStutzki fc Guestertll990() . Using 'syn- 
thesized observations' of the numerically o btained mag- 
netized turbulent clouds, lOssenkopa l)2002t) also showed 
that the X-factor for 12 CO (1-0) linearly increases for 
the column density. They claimed that the X-factor is 
a constant only in a small column density range, and 
this is consistent with our results, although the objects 
(molecular clouds in a galactic disk and a molecular torus 
around the AGN) are very different. 

The line-ratios, such as 12 CO (2-l)/ 12 CO (1-0), are 
often used to infer the physical conditions of molecular 
clouds, mostly in the Milky way (e.g. Oka et al. 1998), 
with the Large- Velocity-G radient (LVG) technique (e.g. 
iScoville fc Solomonlll974j) . We find that the line-ratio is 
not also uniformly observed for the nuclear disk. This 
suggests that the line-ratios observed with large beams 
do not necessarily show the true ratios. In Figure ITU1 the 
maximum line-ratios of 12 CO (l-0)/(2-l), (3-2)/(2-l), (4- 
3)/ (2-1) are plotted for four different beam sizes. One 
finds that the ratios depend on beam sizes, especially 
for higher- J transitions. CO(4-3)/CO(2-l) is different 
by about 20% between results with the finest and lowest 
resolutions. This implies that for a line-ratio obtained 
with a beam size that is comparable to or larger than 
the object (i.e. nuclear disk), the 'true' line-ratio on a 
local scale could be larger than the average by 20%. This 
difference might not be negligible in inferring the physi- 
cal conditions of the ISM using the usual LVG analyses. 
In this sense, we should understand that the currently 
observed molecular line-ratios in external galaxies (e.g. 
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Wall et al. 1993; Zhu et al. 2003) posses some intrinsic 
uncertainty. 

The beam sizes of current radio interferometer observa- 
tions are typically a few arc-seconds, which is comparable 
to or larger than the size of the molecular disk/torus in 
nearby galaxies. Therefore it is impossible to make a 
direct comparison between our models and observations 
on the fine structure of the disk/torus. Here we compare 
our results with radio observations of central regions of 
nearby galaxies on the peak intensity and the line ratio. 
The Virgo high-resolution CO sur vey with the Nobeya ma 
Millimeter- wave Array (NMA) bv lSofue et all l|2001j) re- 
vealed that many normal spiral galaxies have molecu- 
lar gas cores at their central several 100 pc regions, and 
their 12 CO (1-0) intensity is typically 200-400 K km s -1 . 
In Fig. 0Jd, we showed that 12 CO (2-1) intensity is ~ 
3000 K km s _1 on average, and this is larger by one 
order of magnitude than the observed peak intensity of 
normal spiral galaxies Note that, as shown in Fig. Hk>, 
if the nuclear disks/tori have a pc-scale sub-structure, 
the maximum intensity could be much higher than 1000 
K km s _1 . B y combining NMA and the 45m tel escope 
at Nobevama. lSofue et alJ l|2001fl and lKoda et"aD l|2002D 
found a very dense molecular core in the central region of 
NGC 3079, which is Seyfert 2 with a nuclear starburst. 
The 12 CO (1-0) peak intensity of this galaxy is extremely 
high, 3.8 x 10 K km s _1 , and our torus model associ- 
ated with the starburst is consistent with this high peak 
intensity. 

For higher-J transitions, such as CO (3-2) and (4- 
3), the beam sizes of current observations are typically 
~ 20 arc-seconds. Therefore it is hard to make direct 
comparison between the observations and models. Us- 
ing the 15 m James Clerk Max well Telescope (JCMT), 
llsrael fc Baa3p00lll200ll200^ presented J = 2—1, 3—2 
and 4-3 maps of 12 CO , and J = 2 - 1 , 3 - 2 maps of 13 CO 
of 15 nearby spiral galaxies. The peak intensity of 12 CO 
(4-3) distributed fro m 9 K km s" 1 (NGC 278) to 1 019 K 
km s- 1 (NGC 253). IHarrison. Henkel. fc Russell ltl99l 
showed that ico(3— 2) =998 K km s _1 in the central 300 
pc of NGC 253 using JCMT. They also showed that the 
line ratio, 12 CO (3-2)/ 12 CO (2-1)= 1.1 at the center of 
NGC 253. No t e that using the Heinrich- Hertz- Telescope, 
iDumke et all l|200lD showed J C o(3-2) =680 K km s" 1 
and 12 CO (3-2)/ 12 CO (2-1)= 0.7. The line ratio, 12 CO 
(3-2)/ 12 CO (2-1), shown in Fig. Ellis - 1.4, if the beam 
size is comparable to the disk radius. Our results suggest 
that the line ratio in the central 100 pc region is inho- 
mogeneous, and the apparent ratio depends on the beam 
size. To confirm this, we have to await future ultra-high 
resolution observations with the ALMA. 

4. SUMMARY AND DISCUSSION 

In this paper, we made the first 'observational visu- 
alization' of a molecular torus (or disk) around a super- 
massive black hole in an active galactic nucleus under the 
dynamical influence of starbursts through supernova ex- 
plosions. Three-dimensional hydrodynamic simulations 
of the multi-phase gas disk/torus are used as an input for 
the three-dimensional non-LTE radiative transfer calcu- 
lations for 12 CO and 13 CO lines. We found that clumpy 
and turbulent structure of the nuclear disk whose radius 
is about 30 pc can be resolved with a few pc spatial reso- 
lution, which will be achieved with the ALMA for nearby 



galaxies. Our results suggest that the CO-to-H 2 conver- 
sion factor (X-factor) is not uniformly distributed in the 
cental 100 pc, and the X-factor for 12 CO (J =1-0) is not 
constant against the density. On the other hand, 12 CO 
(J =3-2) is nearly constant over the wide range, which 
is ~ 0.3 x 10 20 cm~ 2 (K km s -1 ) -1 , and it can be used to 
estimate molecular gas mass in gas-rich galactic central 
regions associated with active star formation. We also 
found that the observed line-ratios depend on the beam 
sizes, and the average line-ratios (e.g. 12 CO (4-3)/ 12 CO 
(2-1)) do not necessarily represent the 'real' ratios of the 
ISM on a local (parsec) scale. We should be careful to 
use the ratios to infer the physical conditions of the ISM 
if the spatial resolution is not fine enough. 

Although this is the first attempt to investigate the 
three-dimensional structure and dynamics of the ISM in 
the galactic central region, which could be compared di- 
rectly with observations by next-generation radio inter- 
ferometers, our methodology still needs further improve- 
ment in various respects. Firstly, we have not taken into 
account the radiation from the AGN itself. The strong 
UV radiation from the accretion disk could affect both 
the dynamics and excitation of the molecular gases (see 
also Ohsuga & Umemura 2001). We also assumed a uni- 
form UV radiation field. However this uniformity would 
not be the case in reality, because the radiation from 
the massive stars distributed in the clumpy nuclear disk 
should affect the dissociation and excitation of molecules 
as well as the UV heating of the ISM non-uniformly. We 
cannot answer how serious these effects are on the 'syn- 
thesized observations' presented here, before we more 
consistently solve the radiation field. 

The CO abundance, which is assumed constant 
in our calculations, could affect the intensi t y an d 
the CO-to-H2 conversion f actor l)Sakamotol Il996|) . 
lArimoto. Sofue. fc Tsuiimotol (^996) suggested that the 
observed CO intensities depend on the metallicity in 
galaxies, and they found smaller X-factors in the central 
regions than in outer regions of galaxies. In this paper, 
we assumed the solar metallicity, and a smaller X than 
the local value is suggested. Thus the X-factor could be 
even smaller in the nuclear starburst region with higher 
metallicity than the Galactic local value by a factor of 
two or maybe as much as ten. In fact, low values of 
X have been suggested for the central regions of galax- 
ies (e.g. iWall et al.lH993l lReganll2000() . It is certain at 
least that 12 CO (J =1-0) with the standard X-factor in- 
ferred from observations for our Galaxy should not be 
used without caution, especially for the dense ISM with 
strong star forming activities, such as starburst galaxies 
or ultra-luminous infrared galaxies. 

The dense molecular gas can be traced not only by CO 
lines, but also by ionized and atomic carbon, [Cn]158/xm 
and by [Ci] lines. Under some conditions, CO molecules 
are photo-dissociated by UV radiation, but H 2 is still 
present. The neutral carbon may be ionized in the 
photon-dominated region (PDR). Observations with the 
Kuiper Airborne Observatory revealed that [Cn] emis- 
sion is correlated with the 12 CO(1-0) emis sion in gas-rich 
star forming galaxies ijStacev et al.lll99"lj) . The Balloon- 
borne Infrared Carbon Explorer (BICE) also found [Cn] 
emission in the g alactic plane, wh i ch sh ows similar dis- 
tribution to CO l|Nakagawa et alJll998|) . The Antarc- 
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tic Submillimeter Telescope and Remote Observatory 
(AST/RO) has recently revealed that [Ci] emission has 
a spatial extent si milar to low-J CO l ines in the Galac- 
tic central regions (jMartin et al.ll2004j) . This is also the 
case in external galaxies, based on observ ations with the 
Jame s Clerk Maxwell Telescope (JCMT) lftsrael fc Baas] 
12002). These observations sup port the PDR picture 
(e.g. Tielens fe Hollenbachlll985l) in terms of the molec- 
ular clouds illuminated by the nearby OB stars. In the 
present paper, we do not solve the UV radiation field in 
the torus, and photo-chemistry in terms of CO, [Ci], and 
[Cn] is not taken into account. One should note, there- 
fore, that the CO intensities and the conversion factors 
shown in §4 might be modified if the photo-chemistry and 
the UV radiative transfer are fully taken into account. 
Although such calculations are beyond the scope of the 
present paper, we can estimate how important the photo- 
chemical effects are using the published PDR models. 
iMochizuki fc N akagawa (2000) investigated PDR models 
of molecular clouds, solving chemical and thermal equi- 
librium with radiative transfer calculations for simple, 
static molecular cloud s. Their method is bas i cally b ased 
on the PDR model bv lTielens fc Hollenbachl iflQSl . but 
they assume a spherical cloud instead of the plane- 
parallel geometry. The molec ular clouds considered in 
Mochizuki & Nakaeawa (2000) have typically a mass of 
10 - 1O 5 M , a mean hydrogen density 10 2 < n < 10 4 , 
and size 2-20 pc, and these are similar to the clumps in 
our torus models (Wad a fc Normanll2002[) . Mochizuki & 
Nakagawa (2000) found that [Cn] emission originates in 
the outer envelope of the molecular cloud. On the other 
hand, CO emission is concentrated in the inner region of 
the cloud. They also found that luminosity of CO and 
the CO-to-H2 conversion factor are insensitive for the 
far-UV flux over about two orders of magnitude. The 



conversion factor is close to the standard X-factor for a 
relatively massive cloud (> 10Af©). These results imply 
that the UV radiation affects mainly the outer envelopes 
of typical molecular clouds/cores; therefore, the CO in- 
tensity distribution and conversion factor presented in 
the previous section would not be significantly changed 
even if we took into account the photo-chemistry in our 
model. Of course, this must ultimately be verified by 
fully self-consistent radiative transfer calculations with 
chemistry for the clumpy, turbulent medium. 

We showed in this paper that spatial resolution ten- 
times finer than present-day observations in the radio fre- 
quency could drastically change our pictures of the ISM 
in the central regions of the external galaxies. Combin- 
ing numerical techniques like the one presented here with 
the near-future observations by ALMA will be a very 
powerful tool for investigating the unsolved problems as- 
sociated with galactic nuclei, such as the starburst-AGN 
connection and fueling mechanism for the AGNs. 
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Fig. 1. — Surface sections of the density distribution of the gas 
around a SMBH (model A), x-y, y-z, and z-x planes are shown. 
False color is coded according to the local volume density. The 
yellow regions represent n > 10 3 cm -3 . The box size is 64 pc X 
64 pc X 32 pc. 




Fig. 2. — "Tangled network" of the gas around a SMBH with the 
nuclear starburst of Model B. (a) Volume rendering representation 
of the density field. The gray-scale represents 'relative column 
density' or 'relative opacity' distribution, but its absolute value 
does not correspond to any physical scales. The box size is 70 pc 
across, (b) Surface section of the density field of Model B. The 
color is log-scaled volume density (Mq pc - 3 ). 256 3 grid cells are 
used for a (64 pc) 3 computational box. 



8 



Observing the obscuring torus 



-1B0 km/s 





+ 1 80 km/s 













+ 60 kmfe 



Fig. 3. — Channel maps of 12 CO (J = 2 — 1) for a snapshot of 
model A (Figure 1). Six velocity channels from -180 km s _1 to 
+180 km s" 1 arc picked up from 100 channels in ± 500 km s — 1 . 
A viewing angle of the torus 45 degrees from the rotational axis is 
assumed. 




Fig. 4. — Integrated intensity maps of 12 CO (J =2-1) calculated 
from the density, temperature and velocity field of model A shown 
in Figurc[5] Four different 'beam sizes'(shown in white filled circles) 
are assumed to convolve the original map. The unit of axes is the 
number of grid points (i.e. 0.5 pc/grid) 
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Fig. 5. — (a) Volume-weighted phase diagram. The units of 
temperature and density are K and Mq pc~ 3 ~ 40 H cm~ 3 . The 
15 contour levels represent the number of grid-cells in each phase 
(density and temperature), from 10 2 to 10 . (b) Same as (a), but 
for mass-weighted phase diagram. The 15 contour levels represent 
mass in each phase (density and temperature), from 1O 2 M0 to 
10 6 M Q . 
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Fig. 6.— 13 CO integrated intensities of (a) J =1-0, (b) J =3-2, 
and (d) J =4-3 normalized for the intensity of J = 2-1. The input 
model is model B and an edge-on view is assumed. Panel (c) is 
integrated intensity of 13 CO (2-1) (unit is K km s — 1 ). Assumed 
beam-size is shown by a filled circle in each panel. Each box size 
is 64 pc X 64 pc. 




Fig. 7. — Same as Figure but for 12 CO lines. Panel (c) is 
CO (2-1) integrated intensity distribution. 
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Fig. 8. — Two-dimensional distribution of the CO-to-H2 con- 
version factor, X (10 20 cm" 2 [K km s" 1 ]" 1 , for 12 CO (J =1-0). 
Model B with a viewing angle of 45°. Box size is 64 pcx 64 pc. 
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Fig. 9. — CO-to-H2 conversion factor, X [10 20 cm -2 (K km 
s -1 ) - 1 ] as a function of 12 CO integrated intensity (K km s -1 ) 
for four transitions. Each dot represent a grid point of the 2-D 
distribution of the X-factor like FigE] A viewing angle of 45° is 
assumed. 
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Fig. 10. — Dependence of the maximum line-ratios to CO 
(J = 2 — l)on the four "beam sizes", which are 1/8 (open boxes), 
1/4 (open triangles), 1/2 (open diamonds), and 1 (stars) radii of 
the disk. 



